clear all 
clear global

%paths
folderfig = '/Users/wittwer/Dropbox/5_Jason and Milena/data/Matlab_JMP_2R1/main_estimation_submission2R1/3_figures/';
foldermats= '/Users/wittwer/Dropbox/5_Jason and Milena/data/Matlab_JMP_2R1/main_estimation_submission2R1/2_output/';
folderinput = '/Users/wittwer/Dropbox/5_Jason and Milena/data/Matlab_JMP_2R1/main_estimation_submission2R1/1_input';

%Choose model
homedealerL     =1;
dayL            =1;
keepallL        =1;
largetradeL     =0;
tradesizeL      =0;
quantityrobusL  =0;
fejL            =1;

temp=['/MatlabEstimationData_', num2str(dayL), ' _', num2str(homedealerL), ' _', num2str(keepallL),' _', num2str(quantityrobusL),' _', num2str(largetradeL), ' _',num2str(fejL) ];
a_num = importdata(fullfile([folderinput,temp,'.csv']));

data = a_num.data; 
data = sortrows(data, [1,16,17]); %sort Y on id and then dealer 

%structure of data and Y:
    %1)id 
    %2)period id 
    %3)retailer 
    %4)IIROC_SIDE  
    %5)theta 
    %6)large-trade indicator (=1 for trade size > 25 million)
    %7)NaN
    %8)IIROC_YIELD
    %9)venue
    %10)theta 
    %11)theta
    %12)s0j - platform market share for benchmark model, something weird for homedealer model
    %13)rhoj
    %14)xij
    %15)sigma
    %16)bidderid
    %17)countingj (# of line per dealer/id)
    %18)quote on CanDeal jt
    %19) feij
    %20) ctj
    %21) basedealer
    %22) homedealer ID
    %23) base quality (in model extension)
    %24) homedealer quality (in model extension)
    %25) homedealer or not dummy (in model extension)
    %26) countingjH (# of line per dealer-homedealer)
 
    
   
%Fix seed for the random number generator to ensure replicability
rand('twister', 1);

global   eta homedealer quantityrobus day keepall largetrade fej

eta                        = 0;
homedealer          = homedealerL;
quantityrobus       = quantityrobusL;
day                       = dayL;
keepall                 = keepallL; 
largetrade             = largetradeL;
fej                         = fejL;
tradesize               =0;

retailer  = 0;
side = 1;

temp=['/estimation',['_',num2str(side),'_' num2str(retailer),'_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej)]]; 
load(fullfile([foldermats,temp,'.mat'])) 
    
[nanmedian(Mhat); nanmean(Mhat)]
[nanmedian(SE);nanmean(SE)]
nanmedian(p)
    
idxt = 1:size(Y,1);
    
%Back out shocks nu    
theta            = Y(idxt,5); 
price            = Y(idxt,8); 
xij_shocks   = Y(idxt,23) ;

idxotct = idxt; 
n       = size(Y(idxt,:),1);

if side==1
    cutoff = psiSQ + theta(new_period(:,2)) -  xij_shocks(new_period(:,2)) ;
elseif side==-1
    cutoff = psiSQ  - theta(new_period(:,2)) -  xij_shocks(new_period(:,2));
end

if retailer==0
    [f,xi, bw] = ksdensity(price(idxotct),'Bandwidth',2); 
elseif retailer==1
    [f,xi, bw] = ksdensity(price(idxotct),'Bandwidth',1); 
end 

%exclude outliers for the graph
p_1  = quantile(price(idxotct), .001 );
p_99 = quantile(price(idxotct), .999 );

price_ro = price(idxotct);
price_ro = price_ro(price_ro>=p_1 & price_ro<=p_99);

figure
histogram(price_ro, 100,  'FaceAlpha',.1, 'FaceColor', 'black', 'Normalization', 'pdf')
hold on 
plot([nanmean(nanmean(cutoff)) nanmean(nanmean(cutoff))],[0 .275],  'LineWidth',1.5, 'Color', 'black')
xlabel('yield in basis points')
ylabel('probability density')
temp=['/histogram_price_all',['_',num2str(side),'_' num2str(retailer)]]; 
set(gca,'FontSize',19)
set(gcf,'color','w');
saveas(gcf,fullfile([folderfig,temp,'.png']))


%fraction above cutoff:
sum(price(idxotct)>nanmean(nanmean(cutoff)))/size(price(idxotct),1)
  

%end 
